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In  this  work,  the  effect  of  fibre  surface  morphology  on  the  effective  thermal  conductivity  of  the  gas 
diffusion  layer  of  a  polymer  electrolyte  membrane  fuel  cell  is  presented.  Atomic  force  microscopy  was 
used  to  measure  the  fibre  surface  roughness  and  asperity  height  distributions  for  various  fibres  for  Toray 
carbon  paper.  Hertzian  contact  mechanics  was  used  to  determine  individual  micro-contact  areas  and 
thermal  resistances,  and  results  were  compared  with  the  smooth  cylinder  approximation.  The  effective 
thermal  contact  resistance  between  rough  fibres  was  determined  using  resistance  network  theory.  The 
thermal  contact  resistance  and  total  contact  area  were  determined  for  various  angles  of  fibre  orientation 
and  contact  forces;  results  are  presented  as  empirical  formulations.  It  was  found  that  the  effective 
thermal  contact  resistance  is  significantly  affected  by  fibre  roughness  features  when  compared  to  the 
smooth  fibre  case,  which  is  often  used  in  the  literature.  The  analysis  conducted  provides  an  alternative  to 
computationally  expensive  surface  feature  analyses  by  providing  a  tool  which  can  be  used  to  implement 
the  nano-scale  features  of  gas  diffusion  layer  fibres  into  existing  effective  thermal  conductivity  models. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Polymer  Electrolyte  Membrane  (PEM)  fuel  cells  are  electro¬ 
chemical  energy  conversion  devices  that  have  the  potential  to 
generate  electricity  with  zero  local  greenhouse  gas  emissions, 
when  fed  oxygen  from  air  and  hydrogen  gas.  An  understanding  of 
how  thermal  energy  transfer  and  generation  affect  cell  perfor¬ 
mance  and  tunability  is  needed  before  PEM  fuel  cells  can  be  pro¬ 
duced  for  commercial  products.  The  porous  gas  diffusion  layer 
(GDL)  provides  the  main  pathways  for  thermal  conduction  out  of 
the  cell  through  the  solid  phase,  consisting  mainly  of  stacked 
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carbon  fibres.  During  operation,  the  exothermic  electrochemical 
reactions  involved  in  water  generation,  and  as  well  as  ohmic  re¬ 
sistances  caused  by  electrical  current,  form  temperature  gradients 
within  the  cathode  of  the  GDL.  These  thermal  gradients  affect  the 
cell's  relative  humidity,  water  saturation,  and  reaction  kinetics  1  , 
which  are  coupled  to  the  overall  performance  of  the  cell.  It  is 
therefore  crucial  to  understand  how  the  structure  of  the  GDL  can  be 
optimized  to  effectively  transfer  heat  out  of  the  cell. 

Several  thermal  conductivity  models  exist  for  analysing  thermal 
conduction  within  the  PEM  fuel  cell.  Lattice  Boltzmann  methods 
have  been  used  for  determining  the  through-plane  and  in-plane 
thermal  conductivity  of  the  GDL  [2-4].  Wang  et  al.  [2]  used  lat¬ 
tice  Boltzmann  methods  to  analyse  the  effective  thermal  conduc¬ 
tivity  of  unsaturated  GDLs.  They  found  that  the  in-plane  effective 
thermal  conductivity  increases  with  fibre  volume  fraction  and  fibre 
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Fig.  1.  AFM  images  used  in  rough  contact  analysis. 


length,  approaching  a  stable  value  when  the  fibre  length  becomes 
sufficiently  long  [2].  The  authors  also  found  that  the  through-plane 
effective  thermal  conductivity  is  inversely  proportional  to  the 
porosity  of  the  GDL,  which  was  supported  in  lattice  Boltzmann 
studies  presented  by  Yablecki  et  al.  [2,3  .  Yablecki  et  al.  [4]  also 
determined  the  effect  of  GDL  water  saturation  on  the  effective 
thermal  conductivity.  It  was  found  that  the  effective  thermal  con¬ 
ductivity  of  the  GDL  increases  with  increasing  water  saturation, 
with  a  more  significant  effect  in  the  through-plane  direction  rather 
than  the  in-plane  direction  [4].  Yablecki  et  al.  [1]  also  used  an 
analytical  model  to  determine  the  effective  thermal  conductivity  of 
compressed  GDLs.  Their  approach  was  based  on  resistance  network 
modelling,  which  was  first  explored  by  Sadeghi  et  al.  [5,6  .  It  was 
found  that  the  thermal  resistance  in  the  GDL  is  dominated  by  the 
fibre-to-fibre  contacts,  and  the  through-plane  effective  thermal 
conductivity  increases  linearly  with  GDL  compression. 

Although  the  effective  thermal  conductivity  of  the  GDL  has  been 
studied,  it  was  commonly  assumed  that  the  surfaces  of  the  carbon 
fibres  are  smooth.  In  this  study,  the  impact  of  incorporating  realistic 
fibre  surface  morphology,  in  particular  the  circumferential  fibre 
roughness,  in  effective  thermal  conductivity  models  is  investigated. 
The  results  of  this  study  could  be  incorporated  into  other  existing 
effective  thermal  conductivity  models,  further  improving  their  ac¬ 
curacy.  Using  an  analytical  approach  based  on  Greenwood's  rough 
contact  model  [7  ,  the  thermal  contact  resistances  and  contact 
areas  between  rough  fibres  housed  in  the  GDL  are  analysed.  The 
results  obtained  for  various  fibre  orientation  angles  and  contact 
forces  are  compared  to  the  smooth  fibre  cases  found  in  the  litera¬ 
ture  and  are  presented  with  empirical  formulations  representing 
the  mean  and  standard  deviation  of  the  contact  areas  and  thermal 
contact  resistances. 


2.  Experimental  analysis 

Carbon  fibres  housed  in  untreated  Toray  GDL  were  analysed 
using  atomic  force  microscopy  (AFM),  located  at  the  Canadian 
Centre  for  Electron  Microscopy  at  McMaster  University,  Ontario.  A 


Nanoscope  Ilia  Multimode  (Digital  Instruments)  atomic  force  mi¬ 
croscope  was  used  to  image  Toray  TGP-H-120  carbon  fibre  paper 
without  polytetrafluoroethylene  (PTFE)  treatment.  Fibres  from  the 
top  layer  of  the  GDL  were  imaged  to  obtain  the  following  surface 
feature  information:  roughness  in  the  circumferential  and  longi¬ 
tudinal  directions,  surface  area,  and  height  deviations  caused  by 
protruding  asperities  or  irregularities.  Here  it  is  important  to  note 
that  the  longitudinal  roughness  differs  from  the  fibre  waviness.  The 
roughness  in  this  context  can  be  viewed  as  the  deviation  in  the 
distance  from  the  fibre  surface  to  the  central  axis  of  the  fibre  in  the 
direction  considered  (circumferential  or  longitudinal).  The  wavi¬ 
ness  however  can  be  viewed  as  the  path  which  the  central  axis  of 
the  fibre  follows,  where  a  straight  central  axis  would  correspond  to 
a  non-wavy  fibre  [8  . 

Multiple  fibres  were  analysed,  from  which  six  sections  were 
imaged.  These  sections  are  representative  of  the  various  types  of 
surfaces  of  carbon  fibres  exhibited  in  the  AFM  analysis.  AFM  images 
from  two  locations  for  three  fibres  within  untreated  GDL  samples 
are  shown  in  Fig.  1.  The  six  AFM  images  in  Fig.  1  feature  large  and 
small  asperities  (image  b  and  image/,  respectively),  protruding  ir¬ 
regularities  (image  a  and  image  c),  localized  flat-zones  (image  e), 
and  sharp  peaks  (image  d).  A  scanning  frequency  of  1.001  Hz  was 
used  to  obtain  the  images  with  dimensions  of  3  pm  by  3  pm.  The 
image  dimensions  were  determined  by  the  expected  size  of  the 
contact  area  shape  based  on  findings  from  previous  studies 
assuming  smooth  fibre  contact  [1,6];  the  contact  area  does  not 
exceed  this  boundary  for  the  forces  and  orientation  angles 


Table  1 

AFM  image  surface  feature  data. 


AFM  image 

RMS  roughness  [nm] 

Surface  area  pirn2] 

a 

64.794 

10.322 

b 

50.782 

10.427 

c 

95.808 

11.895 

d 

53.864 

10.639 

e 

23.361 

9.807 

f 

75.411 

11.443 
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Fig.  2.  3D  mapping  of  fibre  surface  data  obtained  using  AFM  for  image/. 


considered.  Table  1  includes  statistical  information  regarding  the 
surface  features  for  each  of  the  AFM  images  from  Fig.  1,  including 
the  root-mean  square  (RMS)  roughness  and  surface  area.  Since  the 
rough  fibres  are  nominally  cylindrical  in  shape,  with  an  assumed 
nominal  carbon  fibre  diameter  of  7.32  pm  for  Toray  carbon-paper 
GDLs  [1],  the  roughness  values  in  Fable  1  were  obtained  by 
plane-fitting  the  height  information  to  a  cylindrical  plane  9  .  Since 
the  AFM  images  display  roughness  about  the  cylindrical  surface 
(the  nominal  shape  of  the  fibre),  the  plane-fit  represents  the  height 
deviation  with  respect  to  the  nominal  diameter.  The  plane-fit  was 
performed  in  the  x-direction  of  the  AFM  images  shown  in  Fig.  1  for 
the  roughness  calculation,  as  all  fibres  were  oriented  with  the 
cylinder  axis  parallel  to  the  y-direction. 

Using  the  ASCII  image  data,  a  three-dimensional  mapping  of  the 
fibre  surfaces  were  obtained  using  MATLAB.  Fig.  2  shows  an 
example  of  one  of  the  3D  meshes,  for  AFM  image  /  in  Fig.  1.  The 
image  depicted  in  Fig.  2  includes  the  surface  features  of  a  section  of 
an  individual  carbon  fibre  found  in  the  GDL.  In  the  AFM  images,  the 
x  and  y  axes  represent  the  spatial  co-ordinate  of  the  AFM  image, 
while  the  z  axis  represents  the  height  of  the  surface.  There  is  an 
appreciable  degree  of  roughness  in  the  x  direction  (the  circumfer¬ 
ential  direction  of  the  fibre),  as  opposed  to  the  y  direction  (the 
longitudinal  direction);  this  was  also  observed  in  Ref.  10].  An 
example  fibre  height  profile  for  image  a  is  shown  in  Fig.  3.  Equally 
spaced  slices  were  taken  through  the  length  and  width  of  the  AFM 
images  to  obtain  roughness  profiles. 

While  the  waviness  of  the  fibres  in  the  longitudinal  direction 
may  have  an  important  effect  on  the  nature  of  contact  between  the 
GDL  fibres  and  flat  surfaces,  such  as  the  bipolar  plate  [8  ,  the  focus 
of  this  work  is  to  evaluate  the  fibre-to-fibre  contact  within  the  bulk 
region  of  the  GDL  Nonetheless,  the  fibre  roughness  in  the 
circumferential  direction  was  found  to  be  more  impactful  than  the 
roughness  in  the  longitudinal  direction.  For  example,  Fig.  3(a)  and 


Table  2 

Properties  of  the  GDL  carbon  fibres. 


Average  fibre 

Average  fibre 

Thermal 

Poisson's 

Young's  modulus 

diameter,  d 

length,  1 

conductivity,  kf 

ratio,  v 

of  Elasticity,  Ef 

7.32  |im  [1] 

325  |im  [1] 

120  Witt1  K1 

[1] 

0.3  [14] 

210  GPa  [14] 

(b)  show  the  height  profiles  of  various  cross-sections,  traversing 
through  the  AFM  images  in  the  (a)  circumferential  and  (b)  longi¬ 
tudinal  directions  for  AFM  image  a.  As  can  be  seen  in  Fig.  3(b),  the 
surface  is  relatively  smooth  along  the  fibre  length,  showing  a  single 
asperity  deviation  in  height.  Similar  trends  were  found  for  the 
other  cases.  Throughout  the  various  AFM  images,  it  was  found  that 
the  largest  deviation  in  height  along  the  longitudinal  direction 
caused  by  roughness  (excluding  inclusions  on  the  fibre  surface) 
was  approximately  50  nm,  which  is  negligible  considering  the 
degree  of  roughness  in  the  circumferential  direction  (Fig.  3(a)). 
Height  profiles  for  cross-sections  at  different  positions  along  the 
length  of  the  fibre  are  shown  in  Fig.  3(a).  The  bolded  curve  rep¬ 
resents  the  nominal  diameter  of  the  fibre  [1]  and  was  curve-fitted 
to  each  of  the  respective  AFM  images  using  a  root-mean  square 
approach. 

3.  Methodology 

An  analytical  approach  was  used  to  determine  the  contact  area 
and  thermal  contact  resistance  of  rough  carbon  fibres.  Pairs  of  fibre 
surface  profiles  were  analysed  for  a  range  of  contact  forces  that 
would  be  expected  within  the  GDL  of  an  assembled  fuel  cell.  A 
range  of  fibre  orientation  angles  from  15°  to  90°  was  also  consid¬ 
ered.  The  mechanical  and  thermal  properties  used  in  the  study 
(  ^able  2)  were  assumed  to  be  constant  in  all  directions. 

The  range  of  contact  forces  was  determined  using  the  results  of 
a  previous  study,  where  Yablecki  et  al.  1  developed  analytical 
formulations  for  the  number  of  contact  points  between  adjacent 
layers  of  fibres  in  compressed  GDLs.  These  formulations  were 
derived  using  micro-computed  tomography  data  [11  ,  yielding 
porosity  distributions  for  compressed  GDLs.  Using  the  porosity 
distributions,  this  group  numerically  constructed  GDLs  based  on 
the  nominal  diameter  and  length  of  fibres  housed  in  commercial 
GDLs,  while  assuming  the  fibres  were  preferentially  stacked  and 
allowed  to  overlap  within  adjacent  layers.  The  number  of  contacts 
between  fibres  in  adjacent  layers  was  then  determined  using  the 
numerically  reconstructed  GDLs,  and  presented  as  a  function  of 
orientation  angle  and  average  layer  interface  porosity  [1  .  Using 
the  work  of  [1  ,  we  were  able  to  determine  the  total  number  of 


Fig.  3.  Height  deviation  in  the  circumferential  (a),  and  longitudinal  (b)  directions  for  AFM  image  a. 
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contact  points  (on  average)  between  layers  within  the  GDL  for 
various  GDL  thicknesses.  The  range  of  contact  forces  experienced 
by  the  fibres  in  contact  was  then  determined  using  the  following 
expression: 

Fc  =  (^GDL^GDLWGDL)  j Qayer  0 ) 

where  Fc  is  the  contact  force  experienced  between  touching  fibres 
when  a  GDL  of  length  IqDL,  and  width  wqdl  is  exposed  to  a  through- 
plane  pressure  of  Pqdl-  Equation  (1)  assumes  the  contact  force  is 
evenly  distributed  amongst  the  number  of  contacts  between  the 
two  layers,  defined  as  Qayer.  It  was  found  that  the  maximum  contact 
force  occurs  between  the  layers  with  largest  porosity  values.  The 
largest  contact  force  experienced  between  contacting  fibres  was 
calculated  to  be  0.01  N.  Therefore,  this  value  was  used  as  an  upper- 


bound  when  parameterizing  the  contact  area  and  thermal  contact 
resistance  with  the  contact  force. 

The  focus  of  this  study  is  to  determine  the  impact  of  employing 
fibre  surface  morphology  in  the  determination  of  effective  thermal 
conductivity  of  the  GDL.  To  that  end,  AFM  was  used  to  analyze  the 
fibre  surface  profiles  from  which  an  analytical  methodology  was 
used  to  determine  the  contact  area  and  thermal  contact  resistance. 
Fig.  4  shows  the  flowchart  of  the  algorithm  used  in  this  study  for 
calculating  the  contact  area  and  thermal  contact  resistance  of  GDLs. 

The  first  step  in  the  algorithm  shown  in  Fig.  4  is  the  input  of  the 
fibre  contact  profiles,  which  are  different  than  the  AFM  images 
shown  in  Fig.  1.  The  fibre  contact  profiles  represent  analytical 
mappings  of  each  asperity,  comprised  of  a  diameter  and  position.  It 
was  assumed  that  the  cross-section  of  the  fibres  could  be  arith¬ 
metically  averaged  across  the  3  pm  length  of  each  AFM  image.  The 


Step-wise  compression 


Fig.  4.  Algorithm  used  to  determine  thermal  contact  resistance  and  contact  area  for  rough  contacting  fibres. 
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Fig.  5.  Average  Cross-Sections  with  and  without  curve  fits  for  AFM  images  a,  c,  and  /. 


protruding  asperities  were  curve-fitted  in  each  of  the  average 
cross-sections  with  circular  curves.  The  diameters  and  heights  of 
the  circular  curves  were  obtained  using  a  root-mean  square  curve 
fit,  as  described  in  Ref.  [9  .  The  assumption  that  the  average  cross- 
section  could  be  used  along  the  3  pm  length  of  AFM  images  implies 
that  the  asperities  behave  mechanically  as  small  cylinders.  The 
curve  fits  for  three  of  the  six  AFM  images  from  Fig.  1  are  shown  in 
Fig.  5. 

Once  two  fibre  contact  profiles  have  been  selected  for  analysis, 
the  next  step  is  to  determine  the  number  of  micro-contacts.  A 
micro-contact  is  formed  when  an  asperity  of  one  fibre  comes  into 
contact  with  an  asperity  of  another  fibre.  The  number  of  potential 
micro-contacts  is  the  product  of  the  number  of  asperities  for  the 
two  fibre  contact  profiles  considered.  The  red  asperity  curve  shown 
for  AFM  image  a  is  a  special  case,  which  will  be  described  in  detail 
in  section  4. 

For  this  study,  the  maximum  force  considered  is  0.01  N,  which  is 
not  large  enough  for  the  lowest  asperities  to  come  into  contact  with 
the  tallest  ones.  Therefore,  the  next  step  in  the  analysis  is  to 
determine  the  order  in  which  the  micro-contacts  potentially  come 
into  contact.  For  two  non-parallel  fibres  coming  into  contact  with 
an  angle  of  orientation  ((9),  the  order  of  contact  is  solely  determined 
by  the  height  of  each  individual  asperity.  In  essence,  the  two  as¬ 
perities  with  highest  peaks  will  come  into  contact  first.  If  fibres 
were  to  be  compressed  together  further,  a  second  contact  will  be 
formed  when  one  of  the  two  initially  contacting  asperities  meets 
the  next  tallest  asperity,  and  so  on.  The  distinct  order  of  the  con¬ 
tacts  can  be  described  as  the  ranking  of  the  summed  heights  of  all 
combinations  of  asperities  for  the  two  fibre  contact  profiles 
considered. 

Once  the  contact  problem  has  been  fully  described  using  the 
input  information  of  the  contact  order,  the  angle  of  orientation  is 
considered.  In  this  study,  an  angle  of  orientation  of  90°  means  the 
central  axis  of  the  top  fibre  is  orthogonal  to  the  central  axis  of  the 
bottom  fibre.  The  0°  case  consisting  of  parallel  contacting  fibres  is 
not  considered  in  this  study. 


4.  Contact  mechanics 

The  method  used  in  this  study  for  analysing  the  contact  of  two 
rough  fibres  is  based  on  a  step-wise  compression  algorithm  evolved 
from  Greenwood's  rough  contact  model  [7  .  Since  the  protruding 
asperities  in  this  study  are  fitted  with  circular  curves,  the  asperities 
in  contact  are  treated  as  micro  cylinders,  rather  than  spherical  tips 
as  in  Ref.  [7].  Moreover,  since  there  are  relatively  few  asperities 
within  the  3  pm  by  3  pm  domain,  micro-contact  analysis  is  per¬ 
formed  individually  rather  than  stochastically  via  probabilistic 
functions.  In  this  study,  the  well-established  Hertzian  contact  me¬ 
chanics  formulations  are  used  to  determine  the  micro-contact  area 
and  fibre  normal  displacement  (penetration  displacement).  It  is 
assumed  that  the  individual  asperities  are  smooth,  continuous,  and 
non-conforming  surfaces,  which  may  be  compressed  together 
while  neglecting  friction.  The  Hertzian  equations  are  valid  for  small 
strains,  which  lead  to  stresses  remaining  in  the  elastic  region  for 
these  surfaces.  For  the  Hertzian  analysis  in  this  study,  Poisson's 
ratio  and  elastic  modulus  are  assumed  constant. 

The  three  main  contact  parameters  considered  in  this  study  are 
the  contact  area,  maximum  pressure  on  the  contact  surface,  and  the 
normal  displacement  at  the  center  of  contact.  Due  to  the  curvature 
of  the  contacting  bodies,  the  pressure  within  the  contact  region 
varies.  This  pressure  has  a  maximum  value  where  the  normal 
displacement  is  maximized,  i.e.  the  centroid  of  the  micro-contact 
area,  and  is  a  maximum  when  the  contacting  fibres  are  orthog¬ 
onal.  It  is  important  to  note  that  if  the  maximum  pressure  exceeds 
the  ultimate  strength  of  6.37  GPa  for  carbon  fibre,  the  fibre  will 
fracture  [12  .  Using  a  contact  force  of  0.01  N  does  not  exceed  the 
ultimate  strength  of  carbon,  since  this  force  is  distributed  over  a 
number  of  contacts. 

4  A.  Step-  wise  compression 

The  step-wise  compression  algorithm  begins  with  the  two 
highest  asperities  forming  point  contact.  At  each  step,  the  contact 
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force  is  increased,  and  Hertzian  analysis  is  used  to  determine  the 
contact  area  and  normal  displacement  at  the  locations  of  micro¬ 
contacts.  The  shapes  of  the  micro-contacts  are  generally  elliptical, 
with  the  shape  and  area  being  dependent  on  the  angle  of  orien¬ 
tation  and  the  contacting  cylindrical  asperity  diameters.  As  the 
contact  force  continues  to  increase  incrementally,  a  second  micro¬ 
contact  may  form  when  the  fibre  normal  displacement  is  suffi¬ 
ciently  large.  When  more  than  one  micro-contact  occurs,  the 
increment  in  contact  force  is  distributed  amongst  the  micro¬ 
contacts.  At  the  instance  where  the  second  contact  is  formed,  the 
asperities  forming  the  initial  micro-contact  have  already  been 
compressed  a  certain  amount.  Therefore,  in  order  to  compensate 
for  the  elastic  compression  imposed  up  until  this  point,  the  incre¬ 
ment  of  the  contact  force,  rather  than  the  total  contact  force,  is  split 
up  evenly  amongst  the  micro-contacts  in  the  following  iteration.  In 
this  model,  it  is  assumed  that  the  fibres  are  vertically  compressed 
into  each  other  (without  rotation  about  the  fibre  axis). 

During  each  step,  the  contact  force  is  incremented,  the  algo¬ 
rithm  determines  which  asperity  combinations  form  micro¬ 
contacts,  and  Hertzian  analysis  is  used  to  calculate  the  elliptical 
micro-contact  area  and  normal  displacement  of  the  fibres  for  each 
micro-contact.  The  area  of  each  micro-contact  is  calculated  inde¬ 
pendently,  and  the  effect  of  neighbouring  asperities  is  assumed  to 
be  negligible  in  this  study  since  the  contact  forces  are  small. 
However,  due  to  the  varying  sizes  of  the  asperities,  there  may  be 
virtual  area  overlap  for  asperities  which  are  in  close  proximity  to 
each  other.  Therefore,  the  contact  area  was  examined  to  exclude 
overlapping  micro-contact  regions. 

In  some  cases,  asperity  heights  have  significantly  lower  values 
than  the  majority  of  asperities  ( regular  asperities )  found.  These 
outliers  are  regarded  as  macro-asperities ,  an  example  of  which  is 
shown  in  Fig.  5  (red  curve)  (in  the  web  version).  The  blue  curves 
represent  regular  asperities,  curve-fit  to  each  of  the  protruding 
asperities  in  the  average  cross  sections  for  the  six  AFM  images.  In 
this  model,  macro-asperities  are  utilized  to  handle  fibre  contact 
profiles  which  have  asperities  on  top  of  asperities.  The  goal  of  using 
macro-asperities  is  to  better  approximate  the  contact  area.  For 
example,  in  AFM  image  a  (Fig.  5),  if  the  normal  displacement  were 


to  penetrate  200  nm  into  the  fibre,  the  curvature  of  the  top  regular 
asperity  (curve-fit  with  the  blue  circular  curve)  would  poorly 
describe  the  contact  area  created.  The  model  is  capable  of  deter¬ 
mining  if  the  contact  area  should  be  calculated  using  regular  as¬ 
perities,  or  if  the  regular  asperity  should  be  replaced  with  a  macro¬ 
asperity  to  better  represent  the  contact  scenario.  When  a  macro¬ 
asperity  is  used,  the  initial  compression  of  the  regular  asperity  in¬ 
side  of  the  macro-asperity  will  cause  higher  resistance  to  defor¬ 
mation  than  normal. 

5.  Thermal  analysis 

Once  the  contact  area  is  determined,  the  next  step  is  to  calculate 
the  effective  thermal  contact  resistance  between  the  contacting 
rough  fibres.  The  effective  thermal  contact  resistance  ( TCR )  is 
calculated  using  an  equivalent  resistance  network  representing  the 
constriction  and  spreading  resistances  for  each  individual  micro¬ 
contact  areas.  In  this  study,  it  was  assumed  that  the  only  heat 
transfer  occurring  between  the  fibres  was  steady,  thermal  con¬ 
duction  through  the  carbon  solid  space.  Thermal  conduction  along 
the  length  of  the  fibres  was  assumed  to  be  negligible  in  comparison 
to  the  spreading  and  constriction  resistances.  The  effect  of  PTFE  and 
interstitial  fluids  were  not  considered,  although  the  thermal  con¬ 
duction  is  mainly  through  the  carbon  as  the  thermal  conductivity  is 
much  greater  than  that  of  other  materials  found  in  the  GDL  1  .  Heat 
transfer  caused  by  convection  was  found  to  be  negligible  by 
Ramousse  et  al.  [13]  as  the  velocity  of  fluids  through  the  interstitial 
spaces  and  surrounding  the  fibres  is  miniscule.  Radiative  heat 
transfer  is  negligible  for  temperatures  below  1000  K  2],  which  is 
the  case  for  PEM  fuel  cells  operating  at  temperatures  of  approxi¬ 
mately  80  °C  (353  K).  In  this  study,  the  spreading  and  constriction 
resistance  values  were  assumed  to  be  equivalent,  which  is  in  line 
with  the  studies  conducted  in  Refs.  [5,6,14  . 

Fig.  6  depicts  a  representation  of  the  thermal  and  mechanical 
phenomena  occurring  during  fibre  contact  [  15  .  The  thermal  energy 
is  transferred  between  rough  fibres  through  the  micro-contacts. 
Each  of  the  micro-contacts  have  their  own  respective  spreading 
and  constriction  thermal  resistances  [7  .  The  equivalent  thermal 
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Fig.  6.  Qualitative  representation  of  heat  transfer  between  rough  fibres. 
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resistance  network  used  to  determine  the  TCR  between  two  rough 
fibres  with  temperatures  of  TSink  and  rS0Urce  is  obtained  by  summing 
the  constriction  ( Rco )  and  spreading  ( Rsp )  resistances  in  series,  then 
by  adding  the  sums  ( Rco  +  RSp )  in  parallel  [7]. 

The  formulations  used  to  determine  the  thermal  contact  resis¬ 
tance  for  the  micro-contacts  are  provided  by  Cooper  et  al.  16]. 
Here,  it  is  assumed  that  the  thermal  domain  can  be  treated  as  a 
flux-tube  geometry  considering  the  effect  of  neighbouring  asper¬ 
ities  [15  .  The  formulation  used  is  presented  in  Equation  (2),  where 
fibre  radius  (/fibre)  is  taken  as  3.66  pm  (half  of  the  nominal  fibre 
diameter),  the  thermal  conductivity  of  the  fibre  (kfibre)  is 
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Fig.  7.  Micro-contact  area  versus  contact  force  for  two  AFM  image  a  type  profiles  in 
contact,  for  an  angle  of  orientation  of  (a)  30°,  (b)  45°,  and  (c)  90°. 


Fig.  9.  Total  contact  area  versus  contact  force  for  rough  and  smooth  cases  for  an  angle 
of  orientation  of  45°  (a),  and  total  contact  area  versus  angle  of  orientation  for  rough 
and  smooth  cases  for  a  contact  force  of  0.01  N  (b). 
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Fig.  10.  Average  total  contact  area  contour  plot  versus  angle  of  orientation  and  fibre  contact  force. 
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120  W  m-1  K-1,  and  the  equivalent  micro-contact  radius  (req)  is 
calculated  using  the  micro-contact  area:  req  =  ^/Amc/n.  Once  the 
micro-contact  areas  are  obtained,  the  thermal  contact  resistance 
per  micro-contact  can  be  calculated  using  Equation  (2),  from  which 
the  effective  thermal  contact  resistance  ( TCR )  can  be  computed. 


6.  Results  and  discussion 

The  following  section  includes  the  results  obtained  for  the  an¬ 
alyses  for  contact  area  and  thermal  contact  resistance.  The  algo¬ 
rithm  described  in  the  previous  sections  is  an  analytical  approach 
to  determining  the  contact  area  between  two  rough  GDL  fibres. 
Since  the  surface  morphology  can  vary  between  different  fibres, 
there  is  a  multitude  of  contact  combinations  which  can  occur  even 
with  only  six  fibre  images.  For  example,  the  combinations  studied 
in  this  analysis  could  include  256  cases  formed  by  cross-referencing 
each  of  the  AFM  images,  and  by  vertically  and/or  horizontally 
rotating  fibre  contact  profiles  prior  to  the  analysis.  Each  of  these 
cases  are  analysed  for  various  angles  of  orientation  and  contact 
force.  The  sample  size  for  this  analysis  is  large  enough  to  depict  a 
wide  range  of  fibre  contact  cases  10]. 

6.2.  Contact  area 

Fig.  7  depicts  the  micro-contact  area  versus  contact  force  for  the 
case  of  two  fibre  profiles  with  AFM  image  a  surface  information 


Table  3 

Coefficients  matrix  for  empirical  mean  total  contact  area  and  total  contact  area 
standard  deviation  formulae. 
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coming  into  contact  for  various  angles  of  orientation.  In  these  fig¬ 
ures,  each  of  the  curves  corresponds  to  the  growth  of  an  individual 
micro-contact  area  as  the  contact  force  is  increased.  For  example, 
since  profiles  shown  in  Fig.  7  both  have  two  asperities,  having  them 
come  in  contact  will  lead  to  4  micro-contacts  (four  curves).  In  each 
of  the  three  sub-figures  shown  in  Fig.  7,  some  micro-contacts  are 
formed  after  the  initial  contact.  One  point  of  interest  in  these  fig¬ 
ures  is  the  change  in  slope  of  the  individual  contact  areas.  As  can  be 
seen,  the  rate  of  change  of  the  micro-contact  area  with  contact  force 
decreases  as  new  contacts  are  formed.  The  proportionality  of  the 
micro-contact  area  (Amc)  with  the  contact  force  (Fc)  in  Fig.  7  is  Amc  °c 
Fc //3,  which  is  in  line  with  the  results  of  Archand  [7,17,18J  who  de¬ 
fines  this  proportionality  for  contact  models  which  assume  new 
contacts  can  be  formed  while  existing  contacts  continue  to  grow. 

Fig.  7  also  shows  the  effect  of  altering  the  orientation  angle  on 
the  micro-contact  area;  it  can  be  seen  that  the  values  of  the  micro¬ 
contact  area  decrease  as  the  angle  of  orientation  approaches  90°. 
When  the  orientation  angle  is  90°,  there  is  a  smaller  projection  of 
the  top  fibre  onto  the  bottom  fibre,  resulting  in  smaller  contact 
areas.  When  the  asperities  in  contact  have  equal  diameter  and 
mechanical  properties  (as  in  the  case  depicted  in  Fig.  7(c)),  the 
micro-contact  shape  is  circular.  For  angles  of  orientation  of  30°  and 
45°,  the  micro-contact  area  shape  is  elliptical,  yielding  larger  con¬ 
tact  areas  than  the  orthogonal  case.  The  total  contact  area  is  a 
summation  of  the  individual  micro-contact  areas.  Fig.  8  includes 
total  contact  area  values  for  the  base  cases  depicted  in  Fig.  7,  and  as 
well  as  data  for  two  smooth  fibres  in  contact  (with  no  micro¬ 
asperities).  As  can  be  seen  in  Fig.  8,  the  total  contact  area  for  the 
rough  and  smooth  cases  depicts  a  similar  trend  in  that  the  total 
contact  area  (Atotai)  has  a  proportionality  of  Atotai  00  Fc/3with  the 
contact  force  [7].  The  total  contact  area  for  the  rough  cases  are 
significantly  less  than  the  total  contact  area  for  the  smooth  fibre 
cases,  with  both  cases  having  minimums  at  angles  of  orientation  of 
90°. 

Fig.  9  shows  information  regarding  the  total  contact  area  for 
rough  fibres  in  contact  versus  the  contact  force  at  a  constant  angle 
of  orientation  of  45°  (Fig.  9(a)),  and  versus  the  angle  of  orientation 
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at  a  constant  contact  force  of  0.01  N  (Fig.  9(b)).  In  these  images,  the 
dashed  line  represents  the  smooth  fibre  case;  fibres  with  no  as¬ 
perities.  The  smooth  fibre  approximation  is  typically  considered  in 
existing  effective  thermal  conductivity  analyses  in  the  literature 
[1,5,6,13,14,19,20]  and  provide  a  good  grounds  for  comparison  to 
the  rough  fibre  cases  depicted  in  the  blue  curves.  The  solid,  dark 
curve  represents  the  arithmetically  calculated  average  total  contact 
area  amongst  the  various  fibre  contact  profile  combinations.  The 
shaded  region  represents  the  range  exhibited  throughout  the 
analysis,  defined  by  the  upper  and  lower  bounds  of  the  total  contact 
area  values.  The  total  contact  area  for  the  average  rough  case  and  its 
upper  bound  are  less  than  the  total  contact  area  of  smooth  fibres  in 
contact.  This  is  observed  for  the  range  of  angles  of  orientation  and 
contact  forces.  The  difference  in  the  total  contact  area  for  the  rough 
and  smooth  cases  increases  with  increasing  contact  force,  while  the 
difference  is  minimal  when  the  angle  of  orientation  is  90°.  The 
smooth  fibre  analysis  can  be  repeated  using  Hertzian  mechanics 
and  the  properties  of  the  fibres  provided  in  able  2. 

The  total  contact  area  is  depicted  in  a  contour  graph  as  a  func¬ 
tion  of  angle  of  orientation  and  the  contact  force  in  Fig.  10.  The  data 
displayed  in  Fig.  10  is  the  average  contact  area  amongst  the  various 
fibre  surface  profile  contact  combinations  for  various  angles  of 
orientation  and  contact  force.  As  can  be  seen,  the  contact  area  in¬ 
creases  with  contact  force,  and  decreases  with  angle  of  orientation 
(approaching  the  orthogonal  case  of  90°).  The  change  in  contact 
area  with  angle  of  orientation  is  more  evident  for  larger  contact 
forces.  The  average  measured  rough  contact  area  data,  shown  in 
Fig.  10,  was  surface-fit  to  obtain  an  empirical  formulation  as  a 
function  of  angle  of  orientation  ( 0 )  and  contact  force  (Fc).  The  form 
of  the  empirical  average  contact  area  surface-fit  (Atotai,  avg), 
measured  in  pm2,  is  displayed  in  Equation  (3). 

/  (0,  Fc)  =  ft)  +  Fc  +  @3^  +  @4  OFc  +  fe  F2  (3) 

Since  there  were  many  different  rough  fibre  contact  combina¬ 
tions  considered  in  this  study,  the  standard  deviation  of  the 
measured  contact  area  data  for  each  individual  case  was  also 
calculated  (though  not  included  in  contour  form).  The  form  of  the 
calculated  standard  deviation  in  contact  area  (Atotai,  std)  is  also 
described  in  Equation  (3).  The  constants  for  these  two  different 
functions  are  found  in  Table  3.  The  surface  fits  were  determined 
using  a  least  squares  fit  (MATLAB  2011a),  and  were  fitted  to  the  data 
within  98%  accuracy.  The  surface-fits  for  the  average  contact  area 
and  contact  area  standard  deviation  are  valid  for  the  angle  of 
orientation  range  considered  (15°-90°),  and  contact  force  range  of 
0.0001  N-0.01  N.  Note,  the  units  for  the  input  variables  6  and  Fc  are 
degrees  and  newtons  (N)  respectively.  The  formulations  shown  in 
Equation  (3)  and  Table  3  may  serve  use  to  those  interested  in 
incorporating  fibre  surface  morphology  between  rough  carbon  fi¬ 
bres  within  the  ranges  considered.  Also,  these  formulations  could 
be  used  as  an  input  into  existing  effective  thermal  conductivity 
models  to  incorporate  surface  feature  information  with  respect  to 
the  carbon  fibres.  This  data  would  be  particularly  useful  for  effec¬ 
tive  thermal  conductivity  models  considering  the  effect  of  PTFE, 
binder,  or  interstitial  fluid  at  the  contact  of  rough  fibres  within  the 
GDL.  Users  would  be  able  to  scale  the  amount  of  heat  transfer 
conducting  through  the  solid  carbon,  modifying  their  TCR  to 
include  the  surface  feature  analysis.  Though  the  contact  area  shape 
would  be  needed  to  precisely  determine  the  spreading  and 
constriction  resistances  through  the  micro-contact  areas,  the  total 
contact  area  values  serve  as  a  useful  scaling  parameter  when 
incorporating  roughness  into  existing  effective  thermal  conduc¬ 
tivity  models,  comprising  of  multiple  materials  located  at  the 
contact  region  formed  between  contacting  GDL  fibres. 


6.2.  Thermal  contact  resistance 

The  TCR  data  for  various  fibre  profiles  in  contact  as  a  function  of 
contact  force  and  angle  of  orientation  is  shown  in  Fig.  11(a)  and  (b) 
respectively.  Similar  to  Fig.  9,  the  data  in  Fig.  11  represents  the 
average  rough  TCR  amongst  the  various  fibre  profile  combinations, 
and  as  well  as  the  range  exhibited  through  these  various  combi¬ 
nations,  and  the  analogous  smooth  fibre  contact.  As  shown  in 
Fig.  11(b),  the  TCR  for  the  average  rough  fibres  in  contact  is  signif¬ 
icantly  larger  than  the  smooth  case.  The  range  of  TCR  values  is  large, 
constituting  to  approximately  40%  of  the  average  rough  contact  TCR 
for  contact  forces  greater  than  0.008  N,  increasing  with  decreasing 
contact  force  to  a  value  close  to  90%  of  the  average  rough  contact 
TCR  for  contact  forces  below  0.002  N.  Also,  the  lower  bound  ap¬ 
proaches  the  smooth  fibre  case  for  an  angle  of  orientation  near  90°. 
The  upper  bound  of  the  rough  fibre  TCR  is  defined  by  fibre  profiles 
with  smaller  asperities,  leading  to  smaller  contact  areas,  while  the 
lower  bound  is  defined  by  fibre  profiles  with  larger  asperities  in 
contact.  The  TCR  is  sensitive  to  the  shape  of  the  total  contact  area; 
contact  areas  derived  from  multiple  smaller  micro-contacts  lead  to 
thermal  contact  resistances  which  are  higher  than  those  arising 
from  contact  areas  derived  from  larger  elliptical  shapes. 


Fig.  11.  Effective  thermal  contact  resistance  versus  contact  force  for  rough  and  smooth 
cases  for  an  angle  of  orientation  of  45°  (a),  and  effective  thermal  contact  resistance 
versus  angle  of  orientation  for  rough  and  smooth  cases  for  a  contact  force  of  0.01  N  (b). 


394 


S.J.  Botelho,  A.  Bazylak  /  Journal  of  Power  Sources  269  (2014)  385-395 


Fig.  11(a)  shows  the  effective  thermal  contact  resistance  versus 
the  contact  force.  The  thermal  contact  resistance  for  the  rough  fibre 
case  is  larger  than  the  smooth  fibre  case,  especially  for  lower  con¬ 
tact  forces.  When  the  contact  force  is  significantly  high,  the  amount 
of  compression  of  each  of  the  fibres  into  one  another  becomes 
significantly  large,  and  may  amount  to  the  heights  of  the  asperities 
having  a  minimal  effect  on  the  overall  TCR.  In  essence,  higher 
contact  forces  cause  compression  area  shapes  that  approach  the 
smooth  fibre  case.  The  opposite  is  true  for  lower  contact  forces, 
where  the  surface  roughness  and  asperity  heights  dominate  the 
contact  formed  between  rough  fibres,  leading  to  much  larger  TCR 
values.  This  is  significant  in  the  design  of  GDLs  since  there  will  be 
various  contact  forces  experienced  by  fibres  underneath  the  bipolar 
plates.  In  the  regions  below  the  bipolar  plate  ribs,  there  will  be 
larger  contact  forces,  and  therefore,  the  effect  of  fibre  surface 
morphology  will  be  negligible  with  respect  to  the  computed  TCR. 

A  contour  plot  for  the  thermal  contact  resistance  versus  the 
angle  of  orientation  and  contact  force  is  depicted  in  Fig.  12.  The  data 
in  the  contour  plot  is  calculated  similarly  to  that  of  Fig.  10,  being  the 
arithmetic  average  of  the  various  fibre  profile  contact  combinations 
with  varying  contact  force  and  angle  of  orientation.  The  force  axis  is 
displayed  in  a  logarithmic  fashion,  as  there  are  larger  gradients  in 
thermal  contact  resistance  for  lower  contact  forces.  The  TCR  tends 
to  increase  with  increasing  contact  force  and  as  the  angle  of 
orientation  approaches  the  orthogonal  case.  Also,  the  change  in  TCR 
with  angle  of  orientation  is  more  evident  for  contact  forces  below 
0.001  N;  an  opposite  trend  to  that  shown  in  the  contact  area  con¬ 
tour  in  Fig.  10.  Equation  (4)  shows  the  empirical  form  for  the 
effective  thermal  contact  resistance  average  values  (TCRavg), 
measured  in  K  mW_1,  and  the  standard  deviation  {TCRS tci).  The 
standard  deviation  is  calculated  similarly  to  that  of  contact  area 
standard  deviation,  Atotai,  std-  Empirical  formulations  for  the  TCR 
were  derived  using  MATLAB,  in  a  similar  manner  to  those  derived 
for  the  total  contact  area  found  in  Equation  (3).  The  coefficients  to 
the  two  empirical  surface-fits  are  tabulated  in  fable  4,  and  are  valid 
for  angles  of  orientation  of  15°-90°,  and  contact  forces  of 
0.0001  N-0.01  N.  The  formulations  shown  in  Equation  (4)  and 
Table  4  provide  a  tool  used  to  reproduce  the  analysis  conducted  for 
rough  fibres  in  contact.  It  is  important  to  note  that  the  thermal 


Table  4 

Coefficients  matrix  for  empirically  determined  mean  effective  thermal  contact 
resistance  and  effective  thermal  contact  resistance  standard  deviation. 


/ 

-C0 

«o  x  102 

0L\ 

«2  X  102 

—  0t3  X  101 

a4 

TCRaVg 

6.211 

4.197 

11.9 

6.606 

5.08 

97.42 

TCR  std 

5.307 

-26.120 

10.0 

-2.028 

3.91 

58.43 

contact  resistance  reported  is  solely  due  to  the  conduction  through 
the  solid  space  of  the  contact  carbon  asperities.  Results  from  the 
TCR  study  may  be  validly  used  in  analyses  with  other  materials 
within  the  contact  region  (such  as  PTFE,  binder,  or  interstitial  fluid) 
causing  other  thermal  contact  resistances,  though  care  must  be 
taken  in  summing  the  resistance  values.  The  analysis  conducted 
provides  an  alternative  to  increasing  the  grid  resolution  and  thus 
the  computational  costs  while  still  considering  nano-scale  features 
within  the  fibre  contacts. 

M  Fc)  =  C0  +  +  Fp  +  da*Fp  +  exp(a4Fc)  (4) 

Future  experimental  validation  could  include  the  following  ap¬ 
proaches.  Though  outside  the  scope  of  the  work  presented  here,  the 
fibre-to-fibre  contact  resistance  could  be  measured  directly;  data 
for  which  does  not  currently  exist  in  the  literature.  Also,  this  model 
could  be  incorporated  into  a  macro  model,  such  as  [4,6,21  ,  and 
compared  to  experimentally  measured  bulk  GDL  thermal  contact 
resistances  reported  by  Ref.  22]. 

7.  Conclusion 

In  this  study,  a  method  for  determining  the  contact  area  of  rough 
fibres  in  contact  for  various  angles  of  orientation  and  contact  load 
was  presented.  Fibrous  surface  information  was  obtained  using 
atomic  force  microscopy  for  untreated  Toray  carbon  paper  TGP-H- 
120  GDL  samples  (0%  PTFE  content).  A  modification  of  the  Green¬ 
wood  rough  contact  model,  treating  the  asperities  as  cylindrical 
tips  (rather  than  the  classical  use  of  spherical  surfaces),  was  used  to 
determine  the  contact  area  between  rough  fibres.  It  was  found  that 
the  contact  area  of  rough  fibres  housed  in  the  GDL  changes  with 
angle  of  orientation,  with  minimum  values  when  the  contacting 
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Fig.  12.  Effective  thermal  contact  resistance  versus  angle  of  orientation  and  fibre  contact  force. 
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fibres  are  orthogonal.  The  contact  area  of  rough  fibres  was  found  to 
be  less  than  the  analogous  smooth  fibre  case  for  various  angles  of 
orientation  and  contact  force,  with  larger  deviations  for  higher 
contact  forces.  The  results  were  presented  as  empirical  formula¬ 
tions  which  could  be  implemented  into  existing  effective  thermal 
conductivity  models,  providing  an  alternative  to  computation 
heavy  analyses  considering  nano-scale  features  within  the  GDL 
fibre  contacts. 

It  was  found  that  the  thermal  contact  resistance  decreases  with 
increasing  contact  force,  with  maximum  points  at  angles  of  ori¬ 
entations  of  90°.  The  effect  of  fibre  surface  morphology  on  thermal 
contact  resistance  calculations  is  more  significant  for  lower  contact 
forces,  such  as  those  imposed  on  fibres  underneath  the  bipolar 
plate  channels.  Empirical  formulations  for  the  contact  area  and 
thermal  contact  resistance  were  presented,  and  were  found  to  be 
within  98%  accuracy  of  the  average  contacting  fibre  surface.  The 
goal  of  this  work  was  to  determine  the  impact  of  the  smooth-fibre 
assumption,  and  it  was  found  that  the  circumferential  roughness 
features  do  have  a  significant  impact  on  the  contact  area  and 
thermal  contact  resistance  between  fibres.  The  authors  hope  that 
fuel  cell  modellers  will  now  incorporate  these  roughness  features 
into  their  bulk  (macro)  GDL  models,  which  should  then  result  in 
more  predictive  thermal  conductivity  values  that  can  be  compared 
to  experimental  validation  work. 
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